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We report a systematic study of the weak chemical bond between two benzene molecules. We first 
show that it is possible to obtain a very good description of the C2 dimer and the benzene molecule, 
by using pseudopotentials for the chemically inert Is electrons, and a resonating valence bond wave 
function as a variational ansatz, expanded on a relatively small Gaussian basis set. We employ 
an improved version of the stochastic reconfiguration technique to optimize the many-body wave 
function, which is the starting point for highly accurate simulations based on the lattice regularized 
diffusion Monte Carlo (LRDMC) method. This projection technique provides a rigorous variational 
upper bound for the total energy, even in the presence of pseudopotentials, and allows to improve 
systematically the accuracy of the trial wave function, which already yields a large fraction of the 
dynamical and non-dynamical electron correlation. We show that the energy dispersion of two 
benzene molecules in the parallel displaced geometry is significantly deeper than the face-to-face 
configuration. However, contrary to previous studies based on post Hartree-Fock methods, the 
binding energy remains weak (~ 2kcal/mol) also in this geometry, and its value is in agreement 
with the most accurate and recent experimental findings. [3| 



I. INTRODUCTION 

The intermolecular interaction between benzene rings 
has been a subject of intense theoretical and experimen- 
tal studies in the last two decades[l,[l,0,[l,[l,0|. Indeed 
the intermolecular bonds based on the corresponding tt- 
TT interactions play an important role in many interest- 
ing compounds. For instance, they stabilize the three- 
dimensional structures of biological systems such as pro- 
teins, DNA, and RNA. Moreover, many drugs with a 
specific chemical target utilize these tt-tt interactions and 
the long range forces for their stability. 

In order to understand the mechanism behind those 
attractions, we have considered here the benzene dimer 
as a prototype compound, because both the tt-tt inter- 
actions and the van der Waals (vdW) long range attrac- 
tion are already present and can be studied in a sys- 
tematic way. Despite its simplicity, so far there is no 
general consensus about its equilibrium properties from 
both the theoretical and the experimental side. Indeed, 
it is difficult to determine experimentally the complete 
energy dispersion, and only the total binding energy Dq 
is known, [l], Q with a relatively large experimental er- 
ror due to the weakness of the interactions. On the other 
hand, this compound represents a numerical challenge for 
theoretical methods, because the local density approxi- 
mation (LDA) and other standard treatments based on 
the density functional theory (DFT) are not supposed to 
work well, when dispersive forces are the key ingredient in 
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the chemical bond. Des pite some progress has been made 
recently. (sl. [gI. [qI. [lOl. [Til. Il2l| a general and practical solu- 
tion of this problem is still lacking in the DFT formalism. 
Another family of methods, the accurate post Hartee- 
Fock methods such as CCSD(T), have been extended 
only very recently to a larger basis setQ, since their pro- 
hibitive computational cost has limited their application 
to systems with few electrons and small basis set, and the 
benzene dimer is already at the cutting edge of those ap- 
proaches. As a matter of fact, although the complete ba- 
sis set (CBS) limit can now be estimated more precisely 
in the CCSD(T) framework, the most accurately deter- 
mined binding energy (~ 2.8kcal/mol) of the benzene 
dimer substantially disagrees from the most precise and 
recent measurement [ll] , as also honestly pointed out in 
Ref. 4. Indeed, the CCSD(T) method seems to overbind 
the dimer in the CBS limit. 

Quantum Monte Carlo (QMC) methods are a promis- 
ing alternative to the aforementioned techniques. They 
are able to deal with a highly correlated variational wave 
function, which can explicitly contain all the key ingre- 
dients of the physical system. Their computational cost 
scales favorably with the number of particles N, usu- 
ally as — N^, depending on the method, which makes 
the QMC framework generally faster than the most accu- 
rate post Hartree-Fock (HF) schemes for large enough N . 
Moreover, recent important developments in the QMC 
field allow now to optimize the variational ansatz with 
much more parameters and higher accuracy. In turn 
this can be substantially improved by projection QMC 
methods such as the diffusion Monte Carlo (PMC) [T3| 
and its lattice regularized version (LRDMC). [14 These 
techniques are able in principle to yield the ground state 
energy of the system, since they are based on a direct 
stochastic solution of the Schrodinger equation. How- 
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ever, the well known sign problem affects this kind of 
calculations, and the fixed node (FN) approximation is 
required to make those simulations feasible. Within this 
approximation, it is possible to obtain the lowest vari- 
ational state iPfnIx) of the Hamiltonian with the con- 
straint to have the same signs of a given variational wave 
function ^jq{x). The above condition is applied conve- 
niently in the space representation {x} of configurations 
with given electron positions and spins. It turns out that 
good variational energies can be typically obtained with 
a projection QMC method even starting from a very poor 
variational wave function, the method being clearly ex- 
act in the case when iIiq{x)'iPg{x) > 0,Va;, where V'o(a^) is 
the exact ground state of H. 

Until few years ago, the FN approximation was 
applied[15] to simple variational wave functions obtained 
with basic methods, such as HF or LDA, because for 
large electron systems it was basically impossible to op- 
timize several variational parameters within a statisti- 
cal framework. On the other hand, on small dimer 
systems and even in the single benzene molecule fl7l|. 
it was clearly shown that a highly correlated wave func- 
tion ipQ (x) had to be carefully optimized before applying 
the DMC method with the FN approximation. Other 
examples of the importance of the optimization proce- 
dure have been recently discovered in significant chemi- 
cal systems [l^, showing at the same time that QMC is 
developing quite rapidly and may represent a promising 
tool for future calculations. 

In the present work we report a systematic study 
of the benzene dimer, using the latest developments in 
the QMC framework: an improved optimization algo- 
rithm based on the stochastic reconfiguration (SR), and 
the LRDMC method, which allows to include non 
local potentials (pseudopotentials) in the Hamiltonian 
with a rigorous variational approach. In principle, by 
means of the LRDMC method it is possible to estimate 
Episr = ^'(^"jI,^,/,^^"^ even in presence of pseudopoten- 
tials. Furthermore, a very stable and efficient upper 
bound of Epi^ is obtained by the mixed estimator [19|: 



Elrdmc ~ 
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Elrdmc substantially improves the variational energy 
Eq of the trial wave function ipQ , and is always very close 
to EpN. However, in the case of pseudopotentials, it has 
to be mentioned that EpN obtained with the effective 
Hamiltonian included in the LRDMC is not necessarily 
the lowest variational energy compatible with the signs 
of ipcix). 

The paper is organized as follows: in Sec. [H] we de- 
scribe the variational wave function and its correspond- 
ing basis set. In Sec. lIIIl we introduce the QMC methods 
used. We present some important improvements in the 
SR technique to optimize the energy of a correlated wave 
functions containing several parameters. Moreover, we 
show how it is possible to reduce significantly the lat- 
tice discretization error in the LRDMC method in order 



to improve its efficiency. Finally, in Sec. IIVI we discuss 
the results on the simple but strongly correlated carbon 
dimer, and the more demanding application to compute 
the binding energy of the face-to-face and parallel dis- 
placed configurations in the benzene dimer. 



II. WAVE FUNCTION 

We use the Jastrow correlated antisymmetrized gemi- 
nal power ( JAGP) introduced in Refs. 20 and ,17,, where 
the determinantal part (AGP) is nothing but the par- 
ticle number conserving version of the Bardeen-Cooper- 
Schrieffer (BCS) wave function. The JAGP ansatz is the 
practical representation of the resonating valence bond 
idea, introduced by L. Pauling for chemical systems [2l| . 
and developed also by P. W. Anderson for strongly cor- 
related spin systems [22|. Our variational wave function 
is defined by the product of two terms, namely a Jastrow 
J and an antisymmetric part (^E" — J'^'agp)- The Jas- 
trow term is further split into one-body, two-body and a 
three-body factors (J — J1J2J3) described in the follow- 
ing. All the atomic and molecular cusp conditions are 
fulfilled through the one-body Ji and the two-body J2 
Jastrow factors. The former treats the electron-ion cusp, 
while the latter cures the opposite-spin electron-electron 
cusp. They are both defined by means of a simple func- 
tion u(r) containing only one variational parameter F: 



(2) 



where u'{r) = 1/2 in order to satisfy the cusp condition 
for opposite spin electrons [2^. Then the two-body Jas- 
trow factor reads: 



J2(ri, ...,rAr) = exp 




(3) 



where = jr^— rj| is the distance between two electrons. 
On the other hand the electron-ion cusp condition can be 
satisfied by the one-body term: 



Ji(ri, 



exp 



i,3 



(2Z,f/4^(2Z,)^/^|r.-R,|) 

(4)' 



where Rj arc the atomic positions with corresponding 
atomic number Zj. The reason to take this form for 
the one-body Jastrow factor was inspired by the work of 
Holzmann et a/.jl^] on dense Hydrogen: in the function 
M, the length scaling factor (2Z)^/'' is used to satisfy the 
large distance RPA behavior, whereas the multiplicative 
factor (2^)"^/^ is set by the electron-ion cusp condition: 



< 



dJi 



d\ri - Rj 



>= -ZjJi for |r, - Rj l ^ (5) 



where <> means the angular average. The above relation 
easily follows, since u'{r) = 1/2. 
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Once all the cusp conditions are satisfied, we can 
parametrize the remaining function J3 and the AGP part 
of our resonating valence bond wave function agPi 
and reach the CBS limit for both the full Jastrow factor 
J and the determinantal part, with a Gaussian atomic 
basis set that does not contain any cusp. This rep- 
resents a clear advantage compared with the previous 
parametrization, 17| where it was not even possible to 
satisfy exactly all the electron-ions cusp conditions with 
a finite basis set. Furthermore, this parametrization is 
also particularly useful for interfacing a QMC code with 
standard packages for quantum chemistry calculations, 
which generally use a Gaussian basis set, and therefore 
are not supposed to satisfy any cusp conditions with a 
finite number of basis elements. Obviously this approach 
applies in the same way also for all-electron calculations. 

The AGP geminal function[l3| is expanded over an 
atomic basis set: 



^AGp{r\r 



E 

l.m.a.b 



Al'7</...,(ri)0fc,„.(r^), (6) 



where the indices /, m span different orbitals centered on 
corresponding atoms a, b. In turn, the atomic orbitals 
(j)a.i are expanded with a set of primitive single zeta Gaus- 
sian functions. All the coefhcients and the exponents 
of the gaussians are always consistently optimized. No- 
tice that the largest number of variational parameters are 
contained in the symmetric A matrix, the number of en- 
tries being proportional to the square of the atomic basis 
set size. For this reason, in order to reduce the total num- 
ber of parameters, it is useful to lower the dimension of 
the atomic basis set, by introducing contracted orbitals. 

The three-body J3 Jastrow function takes care of what 
is missing in the one and the two-body Jastrow factors, 
namely the explicit dependence of the electron correla- 
tion on the ionic positions. Therefore, each term in J3 
includes two electrons and one ion interacting each other 
(this is the reason of the name "three-body" ) : 



J3(ri, ...,rAr) = exp ^$,7(rj,rj 



$,7(r„rj) = g,"™V'a,;(ri)V'b,m(rj), (7) 

l,m.a.b 

where the indices I and 771 in the Jastrow geminal <i>j 
indicate different orbitals located around atoms a and 
b, respectively. Again, since all cusp conditions arc al- 
ready satisfied by Ji and J2, in the pairing function 
4>j(ri,rj) we use single zeta gaussian orbitals, i'a.jir) = 
g-zr ^ (simple polynomial in rx, ry, r^.), where fc > is 
an integer and z is the gaussian exponent. The polyno- 
mials are related to the real space representation of the 
spherical harmonics. For instance, to expand J3 up to 
the angular momentum / = 1, we have used two types 
of orbitals, with fc = and fc = 1 respectively. On sim- 
ple dimer compounds we have tested that the inclusion 



of the latter Jastrow orbital is particularly useful for an 
accurate description of the weak vdW interactions. In- 
deed, from a quantum mechanical point of view this type 
of interactions is due to the correlated transition (polar- 
ization) of a couple of electrons from s-wave states local- 
ized around two atoms to corresponding p-wave states. 
Whenever these two atoms are at large distance, we can 
expand J3 for small values of gf'^ , and apply this term 
to a geminal product of two s-wave orbitals. In this way, 
it is clearly possible to describe vdW interactions, pro- 
vided the gaussian basis set used for J3 contains also 
suitable p-wave components. Moreover, we added in the 
J3 pairing function also one body terms, which are the 
product of single zeta gaussian orbitals times a constant 
(i.e. like Qi'^t/ja,!, where c refers to the constant "orbital" 
ipb,c = !)• Thus, our wave function can include a com- 
plete basis set expansion also for the one body Jastrow 
factor. 



III. IMPROVED NUMERICAL METHODS 

In this section we introduce some developments of two 
recently introduced QMC techniques, the SRf^S*! and the 
LRDMCp^ methods, reported in the first and second 
subsection, respectively. The improvements described 
here are of fundamental importance in order to apply 
successfully those methods to realistic electronic systems 
with about 100 valence electrons. 



A. Minimization method 

As described in the previous Section, the JAGP vari- 
ational wave function can contain a large number p of 
non linear parameters {afc}, which are usually difficult 
to optimize for three main reasons, listed below in order 
of difficulty: 

(i) The occurrence of several local minima in the en- 
ergy landscape, leading to the very complex numer- 
ical problem of finding the global minimum energy. 

(ii) The strong dependence between several variational 
parameters. Sometimes, the variation of some 
non linear parameters in the wave function can 
be almost exactly compensated by a correspond- 
ing change of other parameters. This may lead to 
instabilities and/or slow convergence to the mini- 
mum energy. 

(iii) The slow convergence to the minimum energy can 
also be due to simple-minded and/or inefficient it- 
erative methods. 

In the QMC framework, the energy minimization is 
further complicated by the statistical uncertainty, which 
affects all quantities computed, including the optimiza- 
tion target, namely the total energy. Despite these dif- 
ficulties, a lot of progress has been made recently in the 
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energy optimization of highly correlated wave function, 
especially for the alleviation of problems (ii) and (iii) 
[29. [27L [28j. As far as the problem (i) is concerned, the 
solution remains only empirical and relies on the ability 
to find a good starting point of the minimization proce- 
dure. 

In this work we have used a simple improvement of 
the SR method introduced in Ref.^5. for lattice systems, 
and applied later to small atoms 0| and molecules [l7j. 
The SR method has shown to be an efficient and robust 
minimization scheme, although in cases with many vari- 
ational parameters the convergence to the minimum was 
much slower and inefficient for a subset of parameters. 
From this point of view, by using soft-pseudopotentials 
to remove the core electrons, we have experienced a speed 
up in the wave function optimization, because the too 
short wave-length components, responsible of the slow- 
ing down, are no longer present. Moreover, the recent 
methods based on the Hessian matrix [26l [27l. provide 
a further improvement in efficiency, since they allow to 
converge to the minimum energy with fewer iterations. 

Within the SR minimization, the variational parame- 
ters are changed at each iteration: 

a'^. = Uk + Sak 
according to the simple rule: 

Sak^AtJ^Sk'k'h' (8) 

k' 

where At > is small enough to guarantee convergence 
to the minimum, whereas fk = — are the generalized 
forces. The SR matrix s can be any positive definite 
matrix (e.g. if s is the identity matrix one recovers the 
standard steepest descent method), but to accelerate the 
convergence to the minimum and avoid the problem (ii) 
it is much more convenient, as explained in Ref[l3, to use 
the positive definite matrix defined by: 

SkM ^{OkOk')-{Ok){Ok'), (9) 

where the brackets in (C) denote the quantum expecta- 
tion value of a generic operator C over the variational 
wave function tl>G with parameters {ak}- Moreover, Ofc's 
are operators diagonal in the Hilbert space spanned by 
configurations {x}, where electrons have definite posi- 
tions and spins: 

Ok(x) ^ do,,\n\{x\'il:G)\- (10) 

The symmetric matrix s in Eq. [5] has certainly non neg- 
ative eigenvalues because it is just an overlap matrix. In 
the following, we will assume that the matrix s is strictly 
positive definite, as this condition can be easily fulfilled 
by removing from the optimization those variational pa- 
rameters which imply strictly vanishing eigenvalues for 
s. This possibility never occurs in practice, unless the 
wave function has not been efficiently parametrized and 
contains redundant variational parameters. 



At each iteration the various quantities - the ma- 
trix elements s^^k' 1 and the generalized forces - are 
evaluated stochastically over a set of M configurations 
Xi^i = 1,---M, generated by the standard variational 
Monte Carlo method according to the statistical weight 
Tlx — ■'^^[f^- In order to avoid ergodicity problems, 
apparent when the atoms are far apart, we have also in- 
cluded large hopping moves [1^ to the standard Metropo- 
lis transition probability. In the limit M — )■ 00, the sta- 
tistical uncertainty vanishes like 1/\/M, and the above 
minimization strategy certainly converges to some local 
minimum for small enough At and for large enough num- 
ber of iterations. 

In the QMC framework it is obviously important to 
work with a small number A/ of configurations, because 
this number is proportional to the computer time re- 
quired for the optimization. However, though the SR 
method is rather efficient, the statistical noise can dete- 
riorate the stability of the method, especially because the 
matrix s can be ill-conditioned, namely with very small 
eigenvalues, and its inverse can dramatically amplify the 
noise present in the forces. Indeed the SR matrix, even 
when computed with a finite number M of samples, re- 
mains positive definite, but the lowest eigenvalues and 
corresponding eigenvectors can be v ery sensitive to the 
statistical noise. In a previous workpTf . we described a 
simple strategy to work with a well conditioned matrix s, 
by disregarding some variational parameters at each iter- 
ation in the optimization procedure. This method has a 
problem, because sometimes it is necessary to disregard a 
large fraction of the total number p of parameters. More- 
over, we have experienced that removing variational pa- 
rameters from the optimization may be very dangerous, 
as the probability to remain stuck in a local minimum or 
even in a saddle point grows dramatically, especially for 
large p. This occurs even when a relatively small number 
of parameters is not considered in the optimization. 

In order to avoid the above problems, and improve 
the stability of the method, we have modified and sim- 
plified the conditioning of the matrix s. At each step, 
we evaluate the SR matrix with a small bin length 
(M ~ 1000 — 10000), and we regularize it by the sim- 
ple modification of its diagonal elements: 

Sk,k = Sfe,fc(l + e), (11) 

where e can be considered a small Monte Carlo cut-off, 
which can be safely chosen smaller than the average sta- 
tistical accuracy of the diagonal matrix elements Sk,k- In 
this way the modified matrix appears well conditioned 
and without too small eigenvalues. Consequently, the 
improvement in stability can be substantial as shown in 
Fig.dl]) for a simple lattice model test casejl^]. At the 
same time, there is no need to disregard variational pa- 
rameters as in the previous scheme. It is important to 
emphasize that also the modified s matrix is positive def- 
inite, because the sum of two positive definite matrices, 
Si J and £(5, jS, ifsoj. remains a positive definite matrix. 
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As we have already mentioned, this is the only require- 
ment for the iteration in Eq.[5]to converge to a minimum 
ifk = Vfc). Therefore, since all force components fk 
are not biased by the s-matrix modification, by means 
of our approach the exact minimum can be reached for 
arbitrary values of e and Af — )■ oo. 

Obviously other similar regularizations are possible 
and were also adopted elsewhereflE [HI [s^l- For in- 
stance, it is possible to add a simple rescaled identity to 
s {sk,k Sk,k + e), and obtain a well conditioned modi- 
fied matrix with all eigenvalues greater than e. However, 
we have preferred to use the less obvious modification 
in Eq. [Til because in this way the relative change is the 
same for all diagonal elements, which are not deteriorated 
too much in the case they are very small. This is par- 
ticularly useful for the optimization of the present JAGP 
wave function, as it contains some parameters (e.g. the 
Xij in the determinant) ranging in a very tiny interval 
(e.g. within 10"'^ — 10~^) and some others (e.g. the 
exponents Zi in the gaussians) spanning a much wider 
range ( e.g. within 1 — 100). Without the appropriate 
scaling provided by the diagonal elements of the SR ma- 
trix in Eq. [51 an exceedingly small At should be used for 
a stable convergence, which would imply, on the other 
hand, a prohibitively slow convergence. 
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FIG. 1: Optimization of the variational wave function in the 
simple ID Heisenberg model H = J'^^Si ■ Si+i with the 
standard SR (e — 0, open circles), and with the present regu- 
larization (e = 0.001, open triangles). Further details on the 
wave function can be found in Ref. ^7. In the figure, the evo- 
lution of the nearest neighbor spin-spin (n. n. Sz) Jastrow 
parameter is plotted. For each iteration, the forces and the 
SR matrix in Eq. |8] were evaluated over M = 2500 samples, 
whereas At J = 0.125. From this plot it is clear that the SR 
method with e = 0.001 is several order of magnitudes more 
efficient than the standard SR for determining the variational 
parameter with given statistical accuracy. The inset shows 
the first few iterations. 



The present optimization scheme is in practice very 
efficient. For a given bin length M, the SR method be- 
comes optimal for e equal to a finite value, which may 
be even much smaller than the statistical accuracy of the 
matrix elements Sk^k'- In the optimal limit, the statis- 
tical fiuctuations of the variational parameters are sub- 
stantially suppressed without slowing down too much the 
convergence to the energy minimum (see e.g. Fig. [1}. 
Probably, the value e = is optimal only for a noiseless 
infinite precision arithmetic. 

B. LRDMC method with a better o -> limit 

After the energy minimization of a given variational 
wave function ipQ, a substantial improvement in the cor- 
relation energy is obtained by using the DMC method, 
with the so called FN approximation. This method 
allows in principle to determine statistically the low- 
est energy wave function iPfn{x) with the same nodal 
surface as ijjG{x), namely iPfn{x)'iPg{^) ^ (FN con- 
straint). In other words the corresponding energy 

EpN — ^^-^"jl/j^tivT^ minimum possible within 

the FN constraint. Only recently this idea has been 
generalized [3, to include non local potentials in a 
rigorous variational formulation. The LRDMC method 
is based on a lattice discretization of the exact Hamil- 
tonian included in the standard DMC framework. In 
short, the exact Hamiltonian H is replaced by a lat- 
tice regularized one i/", such that — > H for a ^ 0, 
where a is some lattice space which allows to discretize 
the kinetic energy using finite difference schemes, e.g. 

d^^Piy) = V'(y+a)+V'fa-a)-2V-(y) ^ ^i^g^g .g 

trary function. Indeed, our approximate laplacian is: 
A--P /(x,y,z) = 

yy/a^ [p{x + a/2, y, z){f{x + a, y, z) - f{x, y, z)) 
+ p{x - a/2, y, z){f{x - a, y, z) - f{x, y, z))] 
+ X ^ y ^ z. (12) 

where {x,y,z) = r are Cartesian coordinates, and the 
function p is given by 

p(r) = l/(l + Z2|r-R|V4), (13) 

where R is the atom position closest to the electron in 
r, and Z is the largest atomic number considered in the 
system. In particular for the carbon atom, we used Z = A 
throughout this study, as the Is electrons are removed by 
the pseudopotential. The constant rj behaves as l + 0{a^) 
and is introduced to further reduce the error coming from 
the discretization of the kinetic term. 

As pointed out in Refill, an appropriate use of two 
lattice spaces a and a' with fixed irrational ratio a' /a — 
a/ Z^/4: + 1 allows to define in the continuous space 
even for finite a. In the same work, the constant 77 was de- 
termined by requiring that the discretized kinetic energy 
is equal to the continuous one calculated on the state 
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ipG- Here we have found that this requirement is not 
particularly useful for obtaining a very small lattice dis- 
cretization error in the total energy. Indeed as shown in 
Fig. [21 it is much more convenient to define 77 = 1 + Ka^, 
with K determined empirically in order to reduce the sys- 
tematic finite a error. The optimal value oi K — S.2a.u 
(lO.Sa.u.), with a'/a = \/5(\/T0), has been determined 
for the carbon (oxygen) pseudoatom, and can be then 
used also for larger systems containing the same atom, 
as we have done in the forthcoming studies. 

In principle, the LRDMC method allows to calculate 
the expectation value of the Hamiltonian H on the more 
accurate FN wave function ippN- However, this approach 
is rather time consuming because several runs have to be 
performed and some extrapolation is required, which in- 
creases the statistical error by at least a factor of 3. Since 
the LRDMC method - like any other projection method- 
is quite expensive, in the following we have preferred to 
evaluate the simplest upper bound, indicated here by 
Elrdmc, valid for the directly computable mixed aver- 
age Elr.dmc = ^tldV-twT^ ^P^' inequality that 
follows by appl ying the variational theorem on the Hamil- 
tonian i/'' .[lifialal 
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FIG. 2: Energy (Hartree) vs lattice space a for various ways to 
approach the a — > limit. The symbol KEa = KE refers to 
the choice made in Ref. where 77 was obtained by setting 
the lattice regularized kinetic energy equal to the continuous 
one. 



IV. RESULTS 

Before tackling the calculation of the weak interaction 
between two benzene molecules, we studied the effect of 
the basis set on our results, and the size consistency of 
our variational JAGP wave function. The basis set de- 
pendence has been analyzed extremely carefully on the 
carbon and oxygen pseudoatoms, as reported in Sub- 



section IIV Al while the size consistency problem of the 
JAGP ansatz applied to carbon-based compounds is de- 
scribed in Subsection IIVBI We studied the relation be- 
tween the size consistency and the binding energy, com- 
puted for the carbon dimer and the benzene ring. Fi- 
nally, in Subsection IIV CI we report the main results on 
the benzene dimer. 

In all those calculations, we used soft pseudopotentials 
to replace the Is electron pair in the carbon and oxygen 
atoms. The former contains a norm-conserving HF pseu- 
dopotential, generated using the Vanderbilt construc- 
tion [3^, while an ab-initio energy-adjusted HF pseu- 
dopotential is included in the oxygen. In the latter case, 
the effective core potential has been fitted [1^ to repro- 
duce a wide range of HF excitations from the neutral, the 
cation, and the anion atom. The transferability and the 
accuracy of the energy-adjusted pseudopotentials have 
shown to be excellent in a recent systematic study of the 
carbon dimer binding energy [28]. However in this work 
we have not adopted this recent pseudopotentials for the 
Carbon based compound. 



A. VMC/LRDMC basis set dependence on carbon 
and oxygen pseudoatoms 

As shown in Tabs. jH and jlTI, the convergence with 
the basis set appears quite rapid in the carbon and oxy- 
gen atoms. Within the present QMC framework, based 
on the JAGP, there is no need of a large basis set, prob- 
ably because all the cusp conditions can be fulfilled ex- 
actly by the variational wave function, even if it is ex- 
panded over a finite basis set. In this way the polar- 
ization orbitals (e.g. with angular momentum d) have 
not to be included in the QMC ansatz for an accuracy 
of ~ ImH . This is quite remarkable if we consider the 
sensitivity to the basis set commonly observed in conven- 
tional quantum chemistry methods. Indeed, as shown by 
Dunning 37], the contribution of the first polarization d- 
orbital to the correlation energy of the oxygen atom is 
~ 60mH, an effect about two order of magnitudes larger 
than the one reported in Tab. [TTl both for the VMC and 
the LRDMC oxygen atom calculations, where the gain 
in energy (if any) is within the statistical accuracy of the 
simulations (« 0.5mH). In these tables it is interest- 
ing to observe that, while the oxygen is well described 
by a Jastrow-Slater wave function, in the carbon atom 
the AGP plays a crucial role for characterizing the non- 
dynamical correlations, providing an energy gain of about 
lOmH even within the LRDMC method. This shows that 
our approach based on the JAGP is particularly useful for 
generic (saturated and unsaturated) carbon based com- 
pounds. 

In order to extend the calculation to large electronic 
systems, an appropriate contraction of a large primitive 
basis set (up to 6s6p) is important to reduce the dimen- 
sion of the A matrix in the AGP part (Eq. [S]). Notice 
that there is a substantial gain in the LRDMC corre- 



7 



lation energy, by slightly increasing the HF Islp con- 
tracted basis with another contracted s shell. Indeed, 
already the [2s Ip] contraction provides a much better 
LRDMC energy, implying that within our JAGP wave 
function it is possible to improve substantially the nodes 
of the HF Slater determinant, with a little extension of 
the variational freedom. It is important to emphasize 
that we have also optimized the HF determinant in the 
presence of the Jastrow factors described in Sec. [iTl On 
the other hand, as shown in Tab. HI we have obtained 
the HF energy within our general Monte Carlo optimiza- 
tion scheme, even though, in this case, it is obviously 
not necessary to use a statistical method. The LRDMC 
calculation in the HF case was done after optimizing 
the two-body and three-body Jastrow factors, without 
changing the HF determinant. Although the variational 
energy of this Slater-Jastrow wave function is higher than 
the corresponding fully optimized HF -I- Jastrow one, their 
LRDMC energies are the same. This implies that, within 
a single determinant wave function, it is difficult in this 
case to improve the nodes even when the Jastrow and 
the determinantal parts are optimized together. 

Nevertheless, our optimization scheme is very stable 
and reliable and allows to optimize a rather large number 
of variational parameters in a systematic way. Within the 
JAGP ansatz and in particular for the benzene molecule, 
it is extremely important to optimize the wave function 
in order to improve the nodal structure, and obtain a 
good LRDMC total energy. This was previously pointed 
out by two of us, in an all-electron calculation within the 
standard DMC framework. [17] In that case, the DMC 
method provides the same total energies as the LRDMC 
method, used here, since for both methods the FN ap- 
proximation is exactly the same in absence of pseudopo- 
tentials. 



B. Binding energy of C2 and benzene molecule, a 
size consistency study 

Though the [2s Ip] contraction is a rather small basis 
and does not provide the converged result in the total en- 
ergy of the carbon atom, it represents a good compromise 
between accuracy and efficiency, because it can describe 
satisfactorily the chemical bond in all carbon-based com- 
pounds studied, as it is shown in Tab. IIIII 

To this purpose, in this Table we have reported two 
methods to calculate the binding energy. In the standard 
method (method I), we compute the difference between 
the total energy at the equilibrium distance and the sum 
of the energies of the independent fragments for a cho- 
sen atomic basis set. The second method (method II) 
is based on the evaluation of the difference between the 
total energy at the equilibrium distance and the energy 
directly obtained when the constituents of the compound 
are still together, but pulled apart at large distance. 

In the following we show that method II is more appro- 
priate to compute the binding energy within the JAGP 



ansatz. Indeed, the AGP part is the particle conserv- 
ing BCS version only for the total number, not for the 
number in a local sector of the wave function. Therefore, 
if more than one fragment is included in the same AGP 
wavefunction, the number of electrons on each fragment 
is not conserved, and this leads to unphysical charge fluc- 
tuations which are energetically expensive. The Jastrow 
factor can significantly lower the energy, by imposing the 
right occupation number, but the local conservation of 
charge is fully restored only in the CBS limit of the Jas- 
trow expansion. Thus, with a finite basis set in the Jas- 
trow factor, the JAGP wavefunction is clearly more ac- 
curate for a single fragment than for the whole system, 
and method I usually underestimates the binding energy 
of the compound. On the other hand, method II is much 
more accurate, as it includes the cancellation of the finite 
basis set errors in the Jastrow term. 

Moreover, in order to exploit a better cancellation of 
errors, it is important that the energy of the fragments 
at large distance is obtained by iteratively optimizing the 
wavefunction of the compound for larger and larger sep- 
arations of the fragments. In this way, one follows adia- 
batically the fragmentation process, and avoids possible 
spurious energy minima, that may occur in the optimiza- 
tion of a non linear function such as the JAGP. 

For a perfectly size consistent wave function, methods I 
and II should coincide in the CBS limit. The JAGP wave- 
function is perfectly size consistent for fragments which 
are singlet, and with the Jastrow factor in the CBS limit. 
In the case of two singlets at large distance, it is enough 
to define the matrix A of the compound as the sum of 
the two fragments A and B (A = A^"^ + A^'^), with 
an appropriate Jastrow factor freezing the charge in A 
and and B when these two are far apart. In presence 
of unpaired orbitals, e.g. for the triplet Carbon atom, 
size consistency is very difficult to fulfill in general. For 
instance, a singlet S = C2 wavefunction corresponding 
to two entangled Carbon atoms at large distance can be 
obtained only with six independent Slater determinants, 
by appropriately combining the two unpaired p— orbitals 
of each Carbon HF wavefunction, i.e.: 

[S* = 0, A far from B) = 

\p^^ A,Py^ A,P, i B,Py i B) 

+ ^ i APy I APx T B,py t B) 



1 



2V3 
1 

2V3 
1 

2V3 

1 

2V3 



\P^^ A,Py i A,P^^ B,Py i B) 
\P,^ A,Py IA,P, lB,Py^ B) 
\P, iA,Py^ A,P,^ B,Py IB) 

\p, iA,py^ A,p, iB,py^ B) (14) 



where each term in the above expression is a single deter- 
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TABLE L LRDMC ground state energies (Hartree units) for carbon pseudoatom using various basis sets. The LRDMC value 
is a rigorous upper bound for the ground state energy. The hmit a —^ was obtained by using rj = 1 + 3.2a^ in Eq. 1121 
The 2-body Jastrow factor has the form reported in Eq. [3] The Jastrow and the AGP (or HP) geminals are expanded on a 
primitive gaussian basis denoted by (ns mp), where n (m) is the number of s-wave (p-wave) gaussian orbitals. Analogously, 
the number and type of contracted orbitals follow the slash symbol. In particular (6s6p)/[lslp] denotes the standard HF Slater 
determinant. For comparison the HF energy obtained in the CBS limit is —5.319505 Hartree. 



Wave function 


3-body J basis 


AGP basis 


VMC 


LRDMC 


AnP-l-9-hoHv 










AGP+2-body 




(3s3p) 


-5.392 (1) 


-5.416(1) 


AGP+2-body 




(4s4p) 


-5.4066(4) 


-5.4178(3) 


AGP+2-body 




(5s5p) 


-5.4095(3) 


-5.4180(1) 


AGP+2-body 




(6s6p) 


-5.4096(2) 


-5.4181(1) 


AGP+2-body 




(5s5pld) 


-5.4096(2) 


-5.4182(1) 


AGP+2&3-body 


(Islp) 


(5s5p) 


-5.4103(2) 


-5.4181(1) 


HF 




(6s6p)/[lslp] 


-5.3193(3) 


-5.4107(3) 


HF+2&3-body 


(Islp) 


(6s6p)/[lslp] 


-5.3991(3) 


-5.4107(2) 


AGP+2&3-body 


(2s) 


(6s6p)/[2slp] 


-5.4075(2) 


-5.4160(1) 


AGP+2&3-body 


(3s2p) 


(4s5p)/[2s2p] 


-5.4115(1) 


-5.4182(1) 


AGP+2&3-body 


(3s2p) 


(6s6p) 


-5.4113(1) 


-5.4183(1) 



TABLE H: Same as in Tab.® for the oxygen pseudoatom. For a comparison with the reported values, the unrestricted HF, the 
MP2 and the CCSD(T) on the VTZ basis set have total energies of -15.7149, -15.8636, and -15.8822 respectively, calculated 
with Gaussian 03, Revision C.02 '38^. The limit a —> was obtained by using ri = 1 + 10. 8a^ in Eq. 1121 



Wave function 


3-body J basis 


AGP basis 


VMC 


LRDMC 


AGP+2-body 




(2s2p) 


-15.410 (3) 


-15.834(2) 


AGP-f2-body 




(3s3p) 


-15.813(1) 


-15.884(1) 


AGP-f2-body 




(4s4p) 


-15.8611(6) 


-15.8901(3) 


AGP+2-body 




(5s5p) 


-15.8687(4) 


-15.8916(2) 


AGP-f2-body 




(6s6p) 


-15.8685(6) 


-15.8918(3) 


AGP4-2-body 




(5s5pld) 


-15.8679(5) 


-15.8920(3) 


HF+2-body 




(5s5p)/[lslp] 


-15.8674(5) 


-15.8920(3) 



minant, with the orbitals indicated inside the brackets, 
toghcther with the four 2s orbitals (2s '\ A, 2s i A, 
2s] 2s IB). 

The JAGP wavefunction can be perfectly size consis- 
tent even for triplet fragments in the ideal but important 
limit of strong repulsion between electrons in the same 
orbital (strong Hubbard U). In this limit the occupa- 
tion of the same unpaired p— orbital by electrons of op- 
posite spins is forbidden as in the singlet expansion for C2 
(Eq. [T4|) . and the two Carbon Slater determinants, each 
with two unpaired orbitals Px and Py, can be joined in a 
single determinant (AGP singlet), by turning on matrix 
elements such as: 

\Px,Px _ \Py,Py _ \Px,Py _ _\Py,Pm 

■^A,B — '^A,B — ^A,B — ^A,B ' 

where A and B indicate the two Carbon atoms at large 
distance, and these matrix elements are assumed to 
be small compared with the occupied 2s orbitals (e.g. 
A^**'^^ = X^'g — 1). Then, it is simple to show that 
the 6 Slater determinants defining the C2 singlet can be 
obtained with the correct coefficients by a single deter- 
minant AGP wavefunction, provided the double occupa- 
tions of the Px and Py orbitals can be projected out by an 



appropriate Jastrow factor. However, the Jastrow factor 
used here, within the present JAGP expansion, can only 
partially project out double occupation of the same or- 
bital because it depends only on the total electron density 
and not explicitly on the corresponding angular momen- 
tum orbital components. By consequence, the present 
JAGP wavefunction can be only approximately size con- 
sistent in this special case. However, since this loss of 
size consistency is clearly due to a local effect of the cor- 
relation on the same atomic orbital, one expects that 
this contribution should be almost the same both at the 
equilibrium and at large A — B distance, and therefore it 
should affect weakly the chemical bond. 

Within this hypothesis, that looks very well confirmed 
in FigQ, it is possible to obtain a good chemical accu- 
racy (~ O.leV), by using only a single geminal JAGP 
ansatz. Notice that in this picture the LRDMC energy 
appears very smooth and reasonable at large distance, 
even without approaching the energy of two isolated Car- 
bon atoms. This calculation suggests that in the exact 
size consistent framework, which includes many AGP or 
determinants, the total energy should acquire only an ir- 
relevant rigid shift, at least within the LRDMC method. 
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Unfortunately we are not aware of very accurate calcu- 
lation of the full energy dispersion in C2, but the zero 
point energy (ZPE), computed from the data in Fig. ([3]), 
is in very good agreement with the experimental value 
(4.2mi/) 39], clearly supporting the accuracy of our cal- 
culation apart for an irrelevant energy shift. 
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FIG. 3: Energy for two Carbon atoms as a function of their 
distance. We used a (4s5p)/[2s2p] GTO basis set for the 
AGP part, and a (3s2p) uncontracted GTO basis set for the 
Jastrow factor. The atomic basis set convergence has been 
reached within ImH at the VMC level for the molecule at 
the equilibrium distance. The LRDMC and VMC energies 
are not fully size consistent (for large K they should approach 
the energies indicated by the full lines standing below). The 
VMC curve shows a maximum at 7? = 6 and a shallow mini- 
mum at J? = 10, which are almost completely removed by the 
LRDMC energies (within an accuracy of O.leV^). The inset 
shows an expansion of the picture around the equilibrium dis- 
tance. Here a cubic polynomial has been used for fitting the 
data in the range 2.1 < _R < 3. From this interpolation the re- 
sulting VMC [LRDMC] equilibrium distance is 2.357(4) a.u. 
[2.358(5) a.u.], and the ZPE is 4.20(4) mH [4.20(5) mH]. 



The very remarkable outcome of this careful analysis 
is that it is possible to describe well the chemical bond 
in most of the interesting carbon-based compounds, as 
shown in Tab. IIIII As a further confirmation that this 
hypothesis is plausible, we have computed the equilib- 
rium distance of the carbon dimer (Tab. IIV[) using the 
simultaneous optimization of the bond length and the 
variational parameters as described in Ref. By means 
of this technique, based on energy derivatives, we can 
compute the bond length much more accurately than by 
fitting only the total energy around the minimum (e.g. 
in the calculation of the bond length in Fig. [3l a much 
larger statistical error is obtained for this quantity). We 
found a perfect agreement with the experimental bond 
length in the large basis set limit. 



C. The benzene dimer 

As discussed in the previous subsection for two sin- 
glet molecules A and B with electron number Na and 
Nb respectively, the JAGP is size consistent whenever 
the three-body Jastrow factor is optimized in the CBS 
limit. In this way, this term can fully project out the 
charge fluctuations present in the AGP part of the wave 
function, which would erroneously allow a number of elec- 
trons different from N a and Nb even when the molecules 
A and B are at very large distance. In our variational 
wave function the Jastrow geminal (Eq. [7]) is defined 
only on a (3s2p) single zeta gaussian basis set for the 
carbon atom, and a (Is) single zeta for the hydrogen 
atom. Nevertheless, the wave function is very close to be 
size consistent, because the total energy evaluated at a 
fairly large distance, i.e. 12 a.w., is given by Ea+b = 
—75.0825 ± 0.0003i? after a full energy optimization, 
whereas the energy of a single benzene molecule within 
the same basis set is given by Ea = -37.5422 ±0.0002i7, 
i.e. exactly Ea+b/'2 within error bars. Therefore, the 
JAGP ansatz with the chosen basis set is supposed to be 
accurate enough to describe the weak interactions in the 
benzene dimer, as both the basis set convergence and the 
size consistent behavior are taken into account. 

The full dispersion curve of the benzene dimer is re- 
ported in Fig.lO for a face-to-face geometry, together 
with the more accurate LRDMC results. As it is ap- 
parent from this picture, the LRDMC result does not 
change qualitatively the variational outcome, showing 
a very weak dispersion, much less deep if compared to 
the most accurate CCSD(T) results. Our best value 
of the binding energy is 0.5(3) kcal/mol. It is possible 
that the LRDMC method reduces the VMC bond length 
(9 — 10a. M.) by 1 — 2a. u., though an accurate determi- 
nation of this quantity is rather difficult due to the very 
shallow minimum. 

We have extended the calculation to the parallel dis- 
placed geometry (see Fig [61), which has been proposed 
to be the most stable configuration. However, since in 
this case the number of variational parameters is larger 
(~ 10000), we have used partial information of the Hes- 
sian matrix, following the scheme introduced in Ref. (2^ 
to accelerate the convergence of the minimization. In the 
first iteration wc move the parameters along the direction 
gi = f , with s the regularized SR matrix in Eq. 1111 
At this step, the Hessian matrix gives the optimal ampli- 
tude 7i of the parameter change 5d — "figi- Analogously, 
after n iterations the variation of the parameters is given 
by 6d = I]"=i7i5i: where {7ji=i,„ are determined by 
using the Hessian matrix of the last iteration n. After 
changing the variational parameters dk dk + Sdk, a 
new vector gn+i — s~^f is computed with the new wave 
function, and then the procedure is repeated iteratively. 
Notice that a single optimal direction is "collective" in 
the parameter space, as it involves many degrees of free- 
dom. In this way the minimization proceeds in a very 
stable and fast way, as shown in Fig. 21 The main ad- 
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TABLE III: Binding energy (eV) for carbon-based compounds, obtained with the JAGP wave function described in the text 
for a given atomic basis set, reported in the table. The most accurate binding energy is obtained by evaluating the difference 
between the total energy at the equilibrium and at large distance (method II). For the benzene molecule we exploited the 
size consistency of the JAGP ansatz valid for singlet fragments and complete Jastrow factor. Therefore, we considered first 
the fragmentation process CqHq —> 3H2 + 3C2, and then we used the already determined C2 binding energy with method II 
(for II2 the JAGP is clearly size consistent, since it is exact for two electrons). The less accurate method is the standard one 
(method I), obtained by computing the large distance energy by summing the energy of each individual fragment with the same 
basis. In the contracted [2s2pld*], [2s2p*] or [2s Ip*] cases, used generally here in the large system cases, the coefficients of the 
contracted orbitals are assumed to be independent of the angular momentum projection l^. Notice also that the inclusion of 
the polarization d-orbital does not affect the binding of C2 within ImH. The DMC HP binding energy (I) for C2 is 5.66ey |23l. 
The last column refers to the non relativistic value estimated either by experiments or by a very accurate calculation for C2. 



Compound 


3-body J basis 


AGP basis 


# par 


VMC (I) 


LRDMC (I) 


VMC(II) 


LRDMC (II) 


Estimated 


C2 


(3s2p) 


(6s6p)/[2slp*] 


69 


5.806 (16) 


5.946(4) 


6.766(25) 


6.267(4) 


6.36(1) " 


C2 


(3s2p) 


(6s6p)/[2slp] 


74 


5.884(16) 


5.959(4) 


6.862(6) 


6.283(5) 


6.36(1) " 


C2 


(3s2p) 


(4s5p)/[2s2p*] 


95 


5.688(8) 


5.883(4) 


6.910(7) 


6.318(9) 


6.36(1) " 


C2 


(3s2pld) 


(4s5pld)/[2s2pld*] 


136 


5.724(8) 


5.887(4) 


6.893(7) 


6.314(7) 


6.36(1) " 


C2 


(3s2p) 


(6s6p) 


255 


5.763(12) 


5.812(4) 


6.737(5) 


6.289(4) 


6.36(1) " 




(3s2p) 


(6s6p)/[2slp*] 


505 


57.06(3) 


58.105(9) 


59.942(60) 


59.067(8) 


59.24(11) " 



"Ref. 
'■Ref.liO 



TABLE IV: Equilibrium distance of the C2 molecule obtained 
by minimizing the energy of the JAGP with the given basis 
set. The symbols used refer to the ones defined in previous 
tables. 



-150.05 



3-body J basis 


AGP basis 


# par 


R (VMC) 


R (exp) 


(3s2p) 


(6s6p)/[2slp*] 


69 


2.3555(8) 


2.3481 


(3s2p) 


(6s6p)/[2slp] 


74 


2.3559(9) 


2.3481 


(3s2p) 


(6s6p) 


255 


2.3480(6) 


2.3481 



vantage of this method is that the Hessian matrix can 
be calculated in a small basis set, and it is not necessary 
to introduce any further regularization parameter other 
than e IQ-^ (Eq.[ni). 

The optimization of the paraUel displaced benzene 
dimer is rather heavy (about two days on a 64 processor 
SP5 parallel machine) because in every iteration shown 
in Fig. ([5]) a very high statistical accuracy is required, due 
to the so many variational parameters, otherwise all the 
matrices involved in the iteration process (especially the 
large overlap matrix s) are too much noisy. For this rea- 
son we have performed the wave function optimization 
only for two particular geometries reported in Tab. [V] 
From the force components in the two inequivalent di- 
rections, it is clear that the minimum energy occurs at a 
value of i?2 = 7.5±0.2a.u., while Ri is about unchanged. 
The binding energy is 2.2(3) kcal/mol. 



V. CONCLUSION 

In this work, we have devised a QMC framework which 
is able to provide reliable estimates of weak chemical 
bonds, mainly driven by vdW dispersive forces. We used 
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FIG. 4: Total energy (Hartree) of the variational wave func- 
tion during the optimization of all the 10405 variational pa- 
rameters consistent with the chosen basis in the parallel dis- 
placed geometry shown in Fig.(|6}. The case Ri = 7 a.u. and 
R2 = 3.4 a.u. is considered here. In the inset the evolution of 
the variational parameter F (Eq. |3J is shown. 



a highly correlated variational wave function, the JAGP 
ansatz, which contains all the necessary ingredients to 
describe intermolecular interactions: (i) a very high "on 
site" accuracy, through the inclusion of near degeneracy 
correlation effects in the AGP part, (ii) the possibility 
to control the molecular charge distribution through a 
local 3-body Jastrow factor, (lit) the capability to take 
into account the intermolecular correlation, responsible 
for the weak dispersive forces, by means of a "long-range" 
Jastrow term, which connects the molecules (or the frag- 
ments) involved in the interaction. Although the JAGP 



11 



o 

^ 101- 

o 



>> 

O) 

c 

LiJ 



5 - 



- 



'I I I I I I 

V HF 
-A- MP2 

Coupled-Cluster CCSD(T) ext. 
■ VMC JAGP 
• LRDMC+JAGP 



TTJTT- 

6 



TTJTT- 

7 



8 9 10 11 12 13 



R (a.u.) 



FIG. 5: Energy for two face-to-face benzene molecules as a 
function of their distance for different methods. The reference 
was taken at i? = 12. The LRDMC kinetic parameters are 
T] — 1.8, a — 0.5a. It., and a' /a — \/5. The nearest neighbor 
C — C (C — H) distance was set to 2.636 (2.038) a.u. in the 
two molecules. 



ansatz is not size consistent in general, we have shown 
that in the carbon-based compounds analyzed here it is 
possible to obtain accurate and reliable results, by tak- 
ing the calculation of the system in the large distance 
geometry as reference point. 

We have described an improved optimization method, 
based on a proper regularization of the overlap matrix 
in the SR scheme, which can be further boosted by the 
information of the Hessian matrix. With this method it 
is possible to optimize a number of parameters of the or- 
der of 10000. Our fully optimized variational wave func- 
tion has been used as initial guess in projection LRDMC 
calculations. We also found the optimal setting of the 
kinetic parameters in the LRDMC method, in order to 
speed up the diffusion MC simulations with pseudopo- 
tentials. After the optimization step and the LRDMC 
projection, our results arc very weakly dependent on the 
basis set used, at variance with the post-HF quantum 
chemistry methods. 

We studied the face-to-face and displaced parallel ge- 
ometry and energetics of the benzene dimer, which is a 
prototype compound to understand intermolecular dis- 
persive forces. After a full optimization of both the Jas- 
trow and the AGP part, the VMC binding energy re- 
mains in qualitative agreement with the LRDMC result, 
which is supposed to be the most accurate QMC calcu- 
lation. All these findings strongly support the reliability 
of our numerical study. 

The binding of the benzene dimer appears small and 
almost negligible (~ 0.5 kcal/mol) in the face-to-face ge- 
ometry. On the other hand, in the parallel displaced 
configuration where the two molecules are shifted by a 
distance i?i = 3.4a.u., there is a sizable gain in energy. 



TABLE V: Binding energies Ai5 (kcal/mol) and forces 
(kcal/(mol a.u.)) acting on the two independent directions 
Ri and R2 shown in Fig. |6l Energies differences are evaluated 
with respect to the large separation geometry (_Ri = 0,i?2 = 
12a. u.), used also in Fig. O The forces are computed in a 
VMC calculation with the optimized variational wave func- 
tion, and include both Feynman and Pulay contributions. R\ 
and R2 are given in a.u. 



Ri 


R2 


Fi 


F2 


^EvMC 


^Elrdmc 





7 





2.1(2) 


-1.4(4) 


0.2(3) 
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0.1(2) 


0.7(3) 


0.5(3) 


3.4 


7 


0.20(8) 


0.6(1) 


1.4(3) 


2.2(3) 


3.4 


8 


-0.22(6) 


-0.7(1) 


2.0 (3) 


1.8(3) 



which reaches its optimal value of 2.2(3) kcal/mol at 
i?2 — 7 a.u.. Apparently, this is smaller than the most 
recent post HF value ( 2.8 kcal/mol 4]) obtained with 
the CCSD{T) method, after a careful extrapolation to 
the CBS limit. However, by considering the reduction 
of the binding energy due to the zero point vibrational 
energy ZPE {AZPE = 0.37kcal/mol), our result goes 
clearly in the direction of the best experimental estimate 
of the binding energy, which is 1.6 ± 0.2 kcal/mol. [l|. 
The agreement between the experiment and our theoret- 
ical prediction is another striking sign of the capability 
of the QMC techniques to describe accurately not only 
a strong intramolecular bond, but also the very weak in- 
termolecular attractions based on vdW dispersive forces. 





FIG. 6: Geometry of the benzene dimer with the Ri and R2 
distances studied in this work. 
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